Primeros pasos
El primer paso para realizar el análisis es cargar todas las librerías que utilizaremos y posteriormente cargar la base de datos con las observaciones generadas por el sistema de gestión de flota. Adicionalmente se carga un csv que contiene las coordenadas de todos los paraderos de la ruta Metropolitana.
library(tidyverse)
library(sf)
library(lubridate)
library(purrr)
library(kableExtra)
library(leaflet)
library(leaflet.extras)
library(echarts4r)observaciones = read_csv("../tracking.csv")
paraderos = read_csv("../02_Metro_paradas.csv")| event_date | unit_id | route_coloquial_name | lat | lng |
|---|---|---|---|---|
| 19/10/2021 00:19 | 753 | Circuito Metropolitano | 21.02143 | -89.66084 |
| 19/10/2021 00:20 | 753 | Circuito Metropolitano | 21.02143 | -89.66084 |
| 19/10/2021 01:18 | 753 | Circuito Metropolitano | 21.02143 | -89.66084 |
| 19/10/2021 02:17 | 753 | Circuito Metropolitano | 21.02143 | -89.66084 |
| 19/10/2021 02:18 | 753 | Circuito Metropolitano | 21.02143 | -89.66084 |
| Name | title | x | y |
|---|---|---|---|
| 10 | 10 | -89.66343 | 20.96487 |
| 11 | 11 | -89.66343 | 20.96746 |
| 12 | 12 | -89.66344 | 20.97000 |
| 13 SORIANA AV. 128 X AV. MADERO | 13 SORIANA AV. 128 X AV. MADERO | -89.66345 | 20.97234 |
| 14 | 14 | -89.66302 | 20.97529 |
Transformación de clase Sf
Debido a que los data frames anteriores contienen dos columnas para representar la latitud y longitud del punto en cuestión, se puede crear un objeto de tipo sf (simple feature) que permita realizar diversas operaciones espaciales mediante st_as_sf(). Esta agrega una columna adicional al data frame llamada por lo generalmente como geometry que en sí es una lista que contiene las coordenadas del registro. Es necesario declarar el nombre de las columnas en las cuales se encuentran las coordenadas.
observaciones_sf = observaciones %>%
st_as_sf(coords = c("lng", "lat"))
paraderos_sf = paraderos %>%
st_as_sf(coords = c("x", "y"))## Simple feature collection with 609747 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -89.72536 ymin: 0 xmax: 0 ymax: 21.0579
## CRS: NA
## # A tibble: 609,747 x 4
## event_date unit_id route_coloquial_name geometry
## <chr> <dbl> <chr> <POINT>
## 1 19/10/2021 00:19 753 Circuito Metropolitano (-89.66084 21.02143)
## 2 19/10/2021 00:20 753 Circuito Metropolitano (-89.66084 21.02143)
## 3 19/10/2021 01:18 753 Circuito Metropolitano (-89.66084 21.02143)
## 4 19/10/2021 02:17 753 Circuito Metropolitano (-89.66084 21.02143)
## 5 19/10/2021 02:18 753 Circuito Metropolitano (-89.66084 21.02143)
## 6 19/10/2021 02:32 753 Circuito Metropolitano (-89.66084 21.02143)
## 7 19/10/2021 02:43 753 Circuito Metropolitano (-89.66084 21.02143)
## 8 19/10/2021 03:16 753 Circuito Metropolitano (-89.66084 21.02143)
## 9 19/10/2021 04:30 753 Circuito Metropolitano (-89.66084 21.02143)
## 10 19/10/2021 05:29 753 Circuito Metropolitano (-89.66084 21.02143)
## # ... with 609,737 more rows
Un problema que emerge del solo transformar la clase del objeto, es que el CRS (Coordinate Reference System) no se encuentra especificado, lo cual puede llevar a operaciones errones. Para asignar un CRS de manera manual al objeto, se recurre a la función st_set_crs() que recibe un objeto de clase sf y un tipo de proyección basado en el sistema EPSG.
observaciones_sf_crs = observaciones_sf %>%
st_set_crs(value = 4326)
paraderos_sf_crs = paraderos_sf %>%
st_set_crs(value = 4326)## Simple feature collection with 609747 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -89.72536 ymin: 0 xmax: 0 ymax: 21.0579
## Geodetic CRS: WGS 84
## # A tibble: 609,747 x 4
## event_date unit_id route_coloquial_name geometry
## * <chr> <dbl> <chr> <POINT [°]>
## 1 19/10/2021 00:19 753 Circuito Metropolitano (-89.66084 21.02143)
## 2 19/10/2021 00:20 753 Circuito Metropolitano (-89.66084 21.02143)
## 3 19/10/2021 01:18 753 Circuito Metropolitano (-89.66084 21.02143)
## 4 19/10/2021 02:17 753 Circuito Metropolitano (-89.66084 21.02143)
## 5 19/10/2021 02:18 753 Circuito Metropolitano (-89.66084 21.02143)
## 6 19/10/2021 02:32 753 Circuito Metropolitano (-89.66084 21.02143)
## 7 19/10/2021 02:43 753 Circuito Metropolitano (-89.66084 21.02143)
## 8 19/10/2021 03:16 753 Circuito Metropolitano (-89.66084 21.02143)
## 9 19/10/2021 04:30 753 Circuito Metropolitano (-89.66084 21.02143)
## 10 19/10/2021 05:29 753 Circuito Metropolitano (-89.66084 21.02143)
## # ... with 609,737 more rows
Ahora vemos que nuestros objetos tienen asociado el CRS WGS 84.
Creación de buffers
El siguiente paso para nuestro análisis, es crear un buffer para los puntos de los paraderos. La distancia de estos buffers deberá ser de 20 mts. Para esto, se utilizará la función st_buffer() que recibe un objeto de tipo sf y la distancia del buffer.
Sin embargo, ya que que el CRS WGS 84 se basa en grados, es necesario antes realizar una transformación a otro CRS diseñado con base en metros. En ese sentido, se utilizará el EPSG: 6372 Mexico ITRF2008 / LCC utilizando la función st_transform().
paraderos_buffer = paraderos_sf_crs %>%
st_transform(crs = 6372) %>%
select(Name, geometry) %>%
st_buffer(dist = 20)Intersección de observaciones con buffer
Debido a limitaciones del sistema, no es posible conocer el momento exacto en el que una unidad pasó por un paradero. Para solucionar esto, se definió un buffer de 20 mts en cada paradero y se filtrarán únicamente aquellas observaciones que se encuentren dentro de alguno de estos. Para lograr esto se utiliza la función st_filter() que recibe dos objetos de tipo sf y aplica una función de geométrica, en este caso, la intersección st_intersects.
Se debe mencionar que como la proyección del buffer usa el CRS EPSG: 6372, se realizó una reproyección para las observaciones.
observaciones_repro = observaciones_sf_crs %>%
st_transform(crs = 6372)
observaciones_filtered = st_filter(observaciones_repro, paraderos_buffer, .pred = st_intersects)Asociación de observaciones con paraderos
Debido a que existen paraderos cuyos buffers se traslapan entre sí, en muchas ocasiones, una misma observación puede estar contenida dentro de dos o más paraderos, por tanto, se asociará cada observación a el paradero más cercano, evitando así duplicidades en los registros. Lo anterior se realiza con la función st_join() que asocia un paradero a cada observación; la función st_nearest_feature() modifica el comportamiento de la primra función, asegurándose que a cada observación se le asocie al paradero más cercano.
observaciones_nearest = st_join(observaciones_filtered, paraderos_buffer, join = st_nearest_feature)| event_date | unit_id | route_coloquial_name | geometry | Name |
|---|---|---|---|---|
| 19/10/2021 10:13 | 753 | Circuito Metropolitano | POINT (3782833 1047811) | 275 UNIDAD DEPORTIVA KUKULKAN-CBTIS 95 |
| 19/10/2021 10:53 | 753 | Circuito Metropolitano | POINT (3784925 1053110) | 328 |
| 19/10/2021 10:57 | 753 | Circuito Metropolitano | POINT (3784674 1054691) | 335 |
| 19/10/2021 12:15 | 753 | Circuito Metropolitano | POINT (3776638 1053267) | 433 PLAZA DORADA |
| 19/10/2021 13:14 | 753 | Circuito Metropolitano | POINT (3782804 1047744) | 183 CBTIS 95 |
| 19/10/2021 16:26 | 753 | Circuito Metropolitano | POINT (3784112 1048831) | 172 UNIVERSIDAD PEDAGÓGICA |
| 19/10/2021 16:29 | 753 | Circuito Metropolitano | POINT (3784042 1050012) | 291 |
| 19/10/2021 17:24 | 753 | Circuito Metropolitano | POINT (3781877 1056628) | 362 BOXITO KALIA |
| 19/10/2021 18:02 | 753 | Circuito Metropolitano | POINT (3776760 1057372) | 411 GLORIETA LA MESTIZA |
| 19/10/2021 18:24 | 753 | Circuito Metropolitano | POINT (3776225 1053368) | 432 WALMART PENSIONES |
Procesamiento de datos
Una vez que ya se cuenta con una base de datos pre-procesada de las observaciones que se encuentran lo suficientemente cerca de un paradero, el paso siguiente es aplicar un procesamiento más detallado para la final presentación de los hallazgos.
En este punto, los datos sobre las coordenadas de las observaciones ya no son imperativas, de forma que son removidas del análisi para hacer las operaciones más ágiles.
observaciones_pre = observaciones_nearest %>%
tibble() %>%
mutate(
event_date = dmy_hm(event_date)
) %>%
select(-geometry, -route_coloquial_name)Limpieza de datos
Debido a que por cada minuto las unidades envían información acerca de su ubicación, el número de pasajeros, etc., es muy probable que se cuenten con observaciones de una misma unidad que se mantuvo durante más de un minuto en un mismo paradero. Estas observaciones tendrán que ser filtradas del análisis.
Los datos primero se ordenan por el nombre del paradero, luego por la fecha del posteo. Posteriormente, se agrupan de acuerdo a la fecha del posteo, del nombre del paradero y del ID de la unidad. Esto permite observar los momentos en los que la unidad x pasó por el paradero y en el día z.
observaciones_nested = observaciones_pre %>%
arrange(Name, event_date) %>%
group_by(date = as_date(event_date), Name, unit_id) %>%
nest()## # A tibble: 57,409 x 4
## # Groups: unit_id, Name, date [57,409]
## unit_id Name date data
## <dbl> <chr> <date> <list>
## 1 800 1 talleres ado 2021-10-08 <tibble [4 x 1]>
## 2 711 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## 3 743 1 talleres ado 2021-10-08 <tibble [2 x 1]>
## 4 1093 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## 5 731 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## 6 706 1 talleres ado 2021-10-08 <tibble [2 x 1]>
## 7 716 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## 8 726 1 talleres ado 2021-10-08 <tibble [2 x 1]>
## 9 740 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## 10 732 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## # ... with 57,399 more rows
El anterior data frame es un ejemplo de un nesteo, en el cual, una columna de la tabla es una lista de tablas. Con este data frame, se aplicará una función que permite calcular el tiempo que le tomó a la unidad volver a pasar por cada paradero, lo cual finalmente nos permitirá filtrar aquellas observaciones en donde este tiempo sea menor a 5 minutos.
Esta función es una función propia con nombre timediff_na() la cual con base en una columna de fechas, ordenada previamente de más antigua a más reciente calcula el tiempo en minutos entre observaciones consecutivas.
timediff_na = function(df) {
minutos = interval(start = lag(df$event_date, n = 1), end = df$event_date,) %/% minutes(1)
return(minutos)
}Esta función se mapea mediante la función map() de la librería purrr la cual aplica una función determinada a una lista y regresa también una lista.
observaciones_map = observaciones_nested %>%
mutate(
minutos = map(data, timediff_na)
)Posteriormente, se filtran aquellos valores que están en un rango de 1 a 5 minutos, debido a que esto querría decir que a la unidad le tomó entre 1 a 5 minutos dar una vuelta al circuito, lo cual no es factible.
observaciones_filtrado = observaciones_map %>%
unnest(c(data, minutos)) %>%
replace_na(list(minutos = 10)) %>%
filter(!between(minutos, 1, 5))Cálculo de minutos entre unidades por paradero
Una vez depurado el data frame, se procede a realizar un procedimiento análogo, primero nesteando el data frame y luego calculando los minutos entre arribos.
observaciones_final = observaciones_filtrado %>%
ungroup() %>%
select(-minutos, -unit_id) %>%
arrange(Name, event_date) %>%
group_by(date, Name) %>%
nest() %>%
mutate(
minutos = map(data, timediff_na)
) %>%
unnest(c(data, minutos)) %>%
select(event_date, date, Name, minutos)| event_date | date | Name | minutos |
|---|---|---|---|
| 2021-10-08 06:12:00 | 2021-10-08 | 1 talleres ado | NA |
| 2021-10-08 06:29:00 | 2021-10-08 | 1 talleres ado | 17 |
| 2021-10-08 07:31:00 | 2021-10-08 | 1 talleres ado | 62 |
| 2021-10-08 07:55:00 | 2021-10-08 | 1 talleres ado | 24 |
| 2021-10-08 08:14:00 | 2021-10-08 | 1 talleres ado | 19 |
La observación del primer arribo del día por paradero siempre quedará indeterminado, puesto que no hay ningún arribo con el cual compararlo.
Presentación de hallazgos
Con la información ya procesada y debidamente depurada, se estima que para la ruta Metropolitana, durante el 08 de octubre al 22 de octubre, el tiempo de espera promedio fue de 64 minutos.